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The Project for Intercomparison of Land Surface Parameterization Schemes (PILPS) Phase 2(e) 
showed that in cold regions the annual runoff production in Land Surface Schemes (LSSs) is 
closely related to the maximum snow accumulation, which in tutu is controlled in large part by 
winter sublimation. To help further explain the relationship between snow cover, turbulent 
exchanges and runoff production, a simple equivalent model (SEM) was devised to reproduce 
the seasonal and annual fluxes simulated by 1 3 LSSs that participated in PILPS Phase 2(e). The 
design of the SEM relates the annual partitioning of precipitation and energy in the LSSs to 
three primary parameters: snow albedo, effective aerodynamic resistance and evaporation 
efficiency. Isolation of each of the parameters showed that the annual runoff production was 
most sensitive to the aerodynamic resistance. The SEM was somewhat successful in 
reproducing the observed LSS response to a decrease in shortwave radiation and changes in 
wind speed forcings. SEM parameters derived from the reduced shortwave forcings suggested 
that increased winter stability suppressed turbulent heat fluxes over snow. Because winter 
sensible heat fluxes were largely negative, reductions in winter shortwave radiation imply an 
increase in annual average sensible heat. 
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1. Introduction 

The Project for Intercomparison of Land Surface Parameterization Schemes (PILPS) Phase 2(e) 
was designed to evaluate the ability of the current generation of land surface schemes (LSSs) 
used in numerical weather prediction and climate models to represent high-latitude processes. 
Twenty-one LSSs were tested using observed meteorological forcing data for 218 grid cells 
that represented the 58,000 km 2 Tome-Kalix River basin in northern Scandinavia at a spatial 
resolution of 1/4 latitude and longitude (Bowling et al., 2002). . The experimental design and 
summary experiment results are described in Bowling et al. (2002). Nijssen et al. (2002) 
compare the simulated results with observations, and explore the inter-model differences in 
surface temperature, net radiation and runoff production. 

The analyses of Phase 2(e) results presented by Bowling et al. (2002) and Nijssen et al. (2002) 
show that among all the models, annual runoff production is closely related to the maximum 
snow accumulation, which in turn is controlled in large part by winter sublimation. Models 
that detain more of the snow meltwater either as surface water storage or in the soil column, 
tend to generate a more attenuated hydrograph, which better matches observations (Nijssen et 
al., 2002). However, the analyses in Bowling et al. (2002) and Nijssen et al. (2002) are unable 
to reveal conclusive explanations for differences among the models in partitioning of latent and 
sensible heat, surface temperatures, and the effect of stability corrections for turbulent heat 
exchange. In addition, the design of the Phase 2(e) experiments did not tightly control model 
parameters, which complicates isolation of the contribution of specific parameterizations to the 
disparity among simulated results. The intent of this paper is to explore the nature of these 

relationships. 
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Koster and Milly (1997) analyzed a portion of the models from the PILPS Phase 2(a) 
experiment and found that regardless of LSS complexity, they could obtain a reasonable 
estimate of the simulated water balances through the development of simple monthly water 
balance models. The parameters for the monthly models were derived from simulations by the 
individual LSSs. That is, the equivalent models were inferred from simulations by the 
underlying detailed models. These parameters were then used to evaluate the interactions 
between evaporation and runoff processes within and among the LSSs. Gedney et al. (2000) 
used a similar methodology to evaluate predicted changes in surface hydrology as a result of 
climate change. In their experiment, equivalent simple bucket models were inferred for four 
LSSs in three regions, for two different climate scenarios. 

In this paper, we condense the PILPS 2(e) simulation results for each of the models into the 
calibration parameters of a set of simple equivalent models (SEMs), which in turn are used to 
stratify model properties. We expand on the analysis of Koster and Milly (1997) and Gedney 
et al. (2000) by using a simplified monthly model to evaluate the interaction between snow 
cover, evaporation, and runofT processes. Thirteen of the twenty-one models that participated in 
Phase 2(e) returned sufficient model output variables to be used in this analysis. The denvation 
of simplified characterizations of each of these LSS is summarized in Section 3 and Appendix 
A. The results of the SEMs as applied to the Tome-Kalix data sets are presented in Section 4. 
In Section 5, the SEMs are used to explore the results of an inadvertent reduced wind speed and 
shortwave radiation ‘experiment’ conducted as part of PILPS Phase 2(e). 

Five sub-regions of the Tome-Kalix River basins were chosen for analysis, as illustrate in 
Figure 1. The sub-regions are each composed of four 1/4° model gnd cells and were selected to 
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have fairly homogeneous vegetation, as summarized in Table 1. Two of the grid cells are at 
high elevation, and do not contain overstory vegetation; the remaining three grid cells are 
predominantly wooded, with varying annual precipitation and surface water (e.g. lakes and 
wetlands) fractions. Four of the sub-regions correspond approximately with the location of 
gauged sub-basins, two of which were used in other aspects of the PILPS 2(e) experiment 

(Bowling et al. 2002). 

Figure 2 shows the variation of the monthly mean evaporative fraction (ratio of monthly 
evaporation to total available energy, defined in Section 3) with soil wetness, for each of the 
five study regions. Individual models are not identified since the general performance of all 
models is of primary interest. Monthly values were further screened to include only months for 
which the average available energy exceeds an equivalent of 1 mm/day of evaporation and 
snow covered area is less than 50% for each model, typically June through September. The 

unitless soil wetness is defined as: 

|I = (W- W wilt)/ (Wtot -Wwilt) ^ 

where w is soil moisture, w wi „ is the wilting point and w tot is the total soil moisture storage 
capacity. In contrast to our expectations for temperate climates. Figure 2 shows that in the 
majority of the cases the evaporative fraction does not increase with increasing soil water 
availability. Six of the models do show an increase in evaporative fraction with increasing soil 
wetness for at least one of the study regions. However, these six models do not consistently 
show an increase for the same study region. The remaining models demonstrate little change in 
evaporative fraction for increasing soil wetness. There is some tendency for the evaporative 
fraction to decrease with increasing soil wetness due to a slight increase in soil wetness 
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following a decrease in evaporation. This farther confirms that the LSS evaporation is not 
limited by water availability in this low net-radiation environment. 

Figure 2 also shows that the models occupy different equilibrium ranges of soil wetness for 
similar runoff rates, consistent with the findings of Koster and Milly (1997). In general, the 
range of soil wetness values does not change markedly among the five regions, although the 
highest simulated soil wetness values are observed for the region with the greatest annual 
precipitation (Region 2). This is also the region with the highest simulated runoff rates. In 
general, the maximum runoff rates simulated by the models differ more among regions than 
among models, with runoff rates decreasing with decreasing annual precipitation (and 
decreasing annual snow accumulation). Nijssen et al. (2002) have shown that the majority of 
simulated runoff is produced following snowmelt in May and June. Examination of the surface 
and subsurface runoff components returned by each LSS separately (not shown) indicates the 
same relationship with increasing runoff rates with increasing snow precipitation for both 
surface and subsurface flow. This suggests that for the majority of the models the bulk of the 
runoff is produced as quick response surface or subsurface ‘stormflow’ rather than soil 
moisture induced drainage, as is to be expected in a region where runoff is dominated by 
snowmelt, often over frozen soils (soil freezing is represented by all but one of the models). 


3. Simple Equivalent Model Description 

The procedure used closely follows the methods of Koster and Milly (1997) and Gedney et al. 


(2000), with the addition of an energy 


•balance snow model and explicit tracking of the snow- 


covered and snow-free fractions, as illustrated in Figure 3. Tfie SEM, which combines the 
energy-balance snow model and a water balance evaporation and runoff calculation, performs 
energy balance calculations for the snow portion only, and includes three state variables, the 
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root zone soil moisture, total snow water equivalent (SWE), and snow radiative temperature. 
Forcings for the SEM included: a) daily snow and rain precipitation, specific humidity, air 
temperature and downwelling radiation from the Phase 2(e) forcing files, and b) monthly 
interception evaporation, available energy and snow-covered fraction from each LSS. The 

water balance equation is expressed as. 

Aw/ At + AWe/At = P e - E, - R f2 ' 

where w is the total soil moisture, W, is the water equivalent of snow, P, is the effect, ve 

precipitation (rain plus melt minus intercepted evaporation), E, is the total (from snow, 

vegetation and bare ground) evapotranspiration, and R is total runoff. No distinction is made 

between surface runoff and root zone drainage. 

The energy balance of the snow pack is expressed as. 



W e -p.- Cp = 5FF(l-flr)+ (LW-oT;)-Q h -Q le 

(3) 

and 

Q — Pa itfsat (^5 ) _ Qair}/ • 

(3a) 

and 

Q h = P a CPair( T s- T oir)/ ' 

(3b) 


where T s is the snow radiative temperature. The heat capacity of the snow and 
C p and Cp a „ , respectively, the Stephan-Boltzmann constant, ff, as well as the density of air 

and ice ( p, and p , ), are specified constants. Daily average shortwave radiation (SW), longwave 
radiation (L W), specific humidity <?„,>). and air temperature (T„ r „ are input forcings. The snow 
model requires two LSS-dependent parameters, the snow albedo (a) and the effective 
aerodynamic resistance <r„), which are determined by an automatic calibration procedure as 
described below. Snow accumulates when precipitation falls at daily air temperatures less than 
0” C and monthly LSS snow fractions <F S ) exceed 0.01%. The latter restriction is a logistical 
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constraint since the SEM depends on the monthly average snow-covered fraction returned by 
each LSS, the SEM cannot accumulate snow when the LSS predicts a snow fraction near zero. 


. Melt occurs when T s exceeds 0° C. 

Accounting for snow cover, the total evapotranspiration (including transpiration, sublimation 
from the snow pack and bare soil evaporation) is. 


with 


E t = F S *E S + (1-F S )E C 

Ec=Pns*A e +E 0 , 


(4) 


. P — 2 O — ^ sal ^ ' s ^ ^ a,r ^ 

and Cs - 

Here, F s is the monthly snow-covered fraction, E s represents evaporation and sublimation from 
the snow pack and E c includes vegetation transpiration and bare soil evaporation. The effective 
available energy A e is defined by Gedney et al. (2000), as 


Ae = (Rn-G-LE,>a V) 

where R n , G, and LE, are the net radiation, ground heat flux and the latent heat released from 
evaporation of intercepted precipitation, respectively, and X is the latent heat of vaporization. 
The effective available energy was scaled to have the units of a water flux by dividing the 
energy flux by X. The LSS-dependent parameter, p ns , in equation 4, can be viewed as the 
(unitless) model average evaporation efficiency with respect to available energy. The intercept, 
E 0 , is then the average quantity of evaporation which can be sustained by negative sensible heat 

alone. 

One departure from the method of Gedney et al. (2000), is that total evapotranspiration from 
snow-free ground, E c , is solely a function of the total available energy at the ground surface. 
That is, there is no dependence of evaporation on soil moisture. This simplification was 
justified by the preliminary analysis in Section 2, which indicated that evaporation was largely 
independent of soil wetness in the five selected sub-regions. In addition, preliminary 
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experiments with simple models that included a moisture-dependence in the estimation of 
snow-free evaporation yielded a slight decrease in the ability to estimate annual evaporation for 

some models. 

Total runoff is calculated as a piecewise linear function of the root zone soil wetness, in a 
manner similar to the root-zone drainage calculation used by Koster and Milly (1997), as 
follows: 

R(|i)=m R i(ji-pRi) M- <= l l R2; 

R(p)=mR2(|d-pR2) + R(ftR2) ft > ffr.2- 

where |l is the unitless soil wetness as defined by Eq. (1), m R] , m R2 , /to and /to are prescnbed 
constants, and R is runoff in mm/day. Since Eq. (6) is dependent on soil wetness alone, runoff 
is not tied directly to precipitation. This is a reasonable assumption given the seasonal effect of 
snow accumulation in this region, and the lack of runoff response to summer rain events 
(Nijssen et al., 2002). The parameter /to is the soil wetness threshold at which runoff is 
initiated and / i R2 is the threshold for initiation of fast stormflow response. m RI and m R2 are the 
runoff rates (in mm/day) per incremental change in soil wetness. 

Snow-free evaporation and runoff are functions of the updated soil wetness, which makes the 
scheme implicit. The LSS-dependent parameters required for the functions defined by Eqs. (3) 
through (6) (a. r fl , Mr,, M r2 , /to. E 0 , and /L) are derived using objective curve-fitting 
procedures for each model, which are described in Appendix A. A monthly value for the snow 

albedo, is determined for each model for each month. For snow-free evaporation 

calculations, the available energy is based on the net radiation output from each LSS, so a 
snow-free albedo is not needed. Preliminary analysis indicated that the aerodynamic resistance 
must take on different values for stable and unstable conditions (e.g. see Figure A.l). A 
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threshold near-surface temperature gradient is determined for switching between the stable and 
unstable values. The effective aerodynamic resistance, r a , is therefore represented by three 
parameters: r a _stabie, r a _unstabie and AT cr u, as described in Appendix A. 

4. Application of the Simplified Model 

A unique set of equivalent model parameters was estimated for each model, for each LSS for 
each of the five sub-regions, using forcing time series averaged over the four grid cells of each 
time series. The average parameters for all five sub-regions for each LSS are listed m Table 2. 

In order to evaluate the performance of the SEM, the parameter values for each LSS were used 
to integrate Eq. (2) over the ten year time period (1989-98). The ten-year simulation was 
performed twice consecutively to eliminate model initialization effects, with the second ten 
years of simulation used for the analysis. 

The ability of the SEMs to represent the LSSs is summarized m Figure 4 for annual 
evaporation, soil wetness and SWE. The left-hand column shows the two regions with no 
overstory vegetation (Regions 1 and 2); the right-hand column shows the three regions with 
overstory vegetation (Regions 3, 4 and 5). The goodness-of-fit for soil wetness and evaporation 
are presented in terms of mean annual bias. SWE is presented as the cumulative distribution of 
the squared correlation coefficient of LSS and SEM SWE, for each region. For individual 
models, a higher r 2 represents a better fit, but the shape of the CDF yields information about the 
quality of fit for all models in a given region. It can be seen from the mean bias that the simple 
model is able to reproduce the simulated mean annual evaporation and soil wetness fairly well 
for most of the sub-regions. Evaporation appears to be better represented in areas without 
overstory, on average. The SEM fitted to the results of VISA (model M) performed poorly on 
evaporation, particularly in Region 4. VISA frequently had a positive near-surface temperature 
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gradient, which is apparently not well represented by the simple stability scheme. In the 
regions without overstory, ECMWF (model B) was not well represented. In this case, the 
fitting technique tended to overestimate r fl for the stable case, resulting in overestimation of 
sublimation in the spring. Bias in estimation of annual runoff (not shown) is similar to that for 

evaporation. 

Somewhat suprisingly the seasonal cycle of SWE in general was better represented in regions 
with overstory, as shown by the cumulative probability distribution of the (r 2 ) between full 
models and the SEM in Figures 4e and 4f. In this figure the abscissa shows the cumulative 
probability, so that a curve that increases rapidly from zero suggests higher correlation among 
the SEM and full model simulations for a given region. Figure 5 shows the mean seasonal cycle 
of SWE from the LSSs and the simplified model for Regions 2 and 5. The figure shows that 
the simple model is able to capture the differences in seasonal snow both among regions and 
among models, despite errors in the fits. The simplified model both over and underpredicts 
maximum snow extent for different models, primarily due to sensitivity to the stability 

threshold (AT cri ,) in the simplified model. 

5. Characterizing the land surface schemes 
5.1 Analysis of PILPS Phase 2(e) results 

The simple water balance model developed in Section 3 has been shown to replicate simulated 
monthly and annual fluxes fairly well. We can now use the parameters of the simple models to 
investigate the partitioning of energy and precipitation in the LSSs. Because annual 
evaporation is represented without regard for simulated moisture content, the partitioning of 
precipitation in the SEM is controlled by three primary quantities: the effective aerodynamic 
resistance, snow albedo, and evaporation efficiency. The effective aerodynamic resistance is 
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controlled in turn by r a _ slable , r a _ unslable and A T cri , The aerodynamic resistance values found for 
each model, for each region are summarized in Figure 6. There is some tendency for higher 
resistances for the two regions without overstory vegetation, but the difference is not dramatic. 
ISBA (model D) and SEWAB (J) typically have much higher resistances in the stable case. 
These extremes have a limited effect on the differences in model performance, since turbulent 
fluxes are greatly suppressed for all models in the stable case. In the unstable case, CLASS 
(A), HY-SSiB (C), ISBA (D) and SEWAB (J) tend to have the highest resistances, with direct 
consequences for total snow accumulation, as shown in Figure 7. In Figure 7a, average annual 
total melt (a surrogate for maximum accumulation) is plotted versus the average effective 
aerodynamic resistance. Regions 2 and 5 are represented by different colors. Because Figure 
7a shows the average r fl , it incorporates the effect of both r fl values, as well as A T crit .. Although 
CLASS does not have the highest resistances for Region 2, its effective resistance is the highest 

due to the relatively low value of A7 1 cr n ("3.2 C). 

Figure 7b shows a similar figure for average snow albedo versus average annual total melt. 
Although the relationship between albedo and total snow accumulation is less well defined, 
there is still a noticeable effect in some cases. For example, NOAH (model H) and VISA (M) 
have similar effective resistances in Region 5, but large differences in melt. This difference is 
explained in part by the large difference in average albedo between the two models. 

The influence of the evaporation efficiency is explored in Figure 8. This figure shows the 
mean annual evaporation efficiency versus annual available energy for Regions 2 and 5 for 
months with snow cover less than 10%. The vectors show the direction and magnitude of 
interannual variability (as represented by the standard deviation) of each of the quantities. The 
angle therefore indicates the relative influence of each quantity on the interannual variability of 
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summer evaporation from each model. For Region 2, these angles vary between 10° and 27°. 
For Region 5, the angles vary between 16° and 43°. This indicates that in all cases, the 
interannual variability in annual evaporation is explained in greater proportion by the changes 
in atmospheric forcings, rather than by differences in model parameterization. 

Roster and Milly (1997) discuss model differences in terms of the average evaporation 
efficiency and the fraction of active soil moisture over which runoff occurs. They found that 
these quantities explain 94 percent of the intermodel variability in evaporation, thus clarifying 
that different evaporation rates in the LSSs resulted from different formulations of canopy 
interception and differences in the relative shape and position of the effective runoff and 
evaporation functions. For the Phase 2(e) results, our analysis indicates that in cold regions 
available energy rather than the soil moisture range controls intermodel differences in summer 
evaporation. As demonstrated by Nijssen et al. (2001), the intermodel difference in net 
radiation over the summer is equally attributable to differences in surface radiative temperature 

and land surface albedo. 

To further explore the sensitivity of mean annual runoff to the aerodynamic resistance, albedo, 
snow cover fraction and evaporation efficiency, a sensitivity analysis was conducted in which 
each of these parameters was systematically fixed to the same value for all models for all years. 
In the case of snow cover fraction, the fixed value(s) were the average over all models of the 
monthly snow cover fractions. Snow albedo likewise was fixed at the average over all models, 
but was not changed seasonally. Figure 9a shows annual ninoff for Region 2 predicted by each 
of the SEMs when snow cover fraction, albedo, resistances and evaporation parameters were 
fixed for all models and regions, listed in Table 3 (the ‘fixed 1 case). There are still some 
differences in simulated runoff among models due to differences in runoff foimulation. 
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However, these differences are small compared to the differences between the fixed case and 
the original runoff estimation using the simplified model (Figure 9b). Figures 9c - 9f show the 
percent difference between the perturbed run and the fixed run for the following perturbations: 
(1) fixed R .’ s and AT„,„ (2) fixed E 0 , and /L, (3) fixed albedo and (4) fixed snow cover 
fraction. The percent difference is calculated as the difference between the perturbed case and 
the original simulation, normalized by the difference between the fixed case and the original 
simulation and multiplied by 100%. For Region 2, the difference between the fixed and 
original simulations is largely explained by the tnter-model differences in effective 
aerodynamic resistance, and to a lesser extent by differences in evaporation efficiency. As 
shown in Figure 10, the situation is less clear for Region 5, which has overstory vegetation. In 
Region 5, aerodynamic resistance and evaporation efficiency still show the most sensitivity for 
most models, but the results are less consistent, and the differences in runoff formulations seem 

to play a larger role. 

5.2 Sensitivity Analysis: results of the Phase 2(e) base run experiment 

As indicated in Bowling et al. (2001), the original forcing data distributed to the participants 
had a low bias in shortwave radiation. The bias resulted from use of a general 
parameterization, based on the method of TV A (1972), to reduce radiation at the top of the 
atmosphere to clear sky radiation at the ground surface. A re-calculation of the shortwave 
forcing using the method of Eagleson (1970) with the turbidity factor estimated from station 
data in the P1LPS 2(e) region, resulted in an average annual increase in shortwave radiation 
basin-wide of 36 percent. The revised data were distributed for a re-run experiment. The re- 
nin data also included a corrected gridded wind field and hourly disaggregated precipitation. 
Correction of the error in the gndding of the wind speed dataset resulted in an increase in wind 
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speed of about 10 percent on average, with greater spatial variability than the shortwave 
radiation correction (in fact, although the domain average increased, some parts of the domain 
experienced an average decrease in wind speed). As discussed in Bowling et al. (2002), the 
disaggregation to hourly precipitation did not affect most of the models substantially. 

Twenty of the LSSs involved in the Phase 2(e) experiment returned results for both the re-runs 
(higher shortwave) and the base runs (lower shortwave), thus leading to an interesting, if 
unplanned, sensitivity experiment. For simplicity of further discussion these are refered to as 
experiment (1) and experiment (2), respectively. All of the analysis presented so far has 
utilized the results of experiment (1). The primary result of the change between experiment (1) 
and (2) was a decrease in annual sensible heat flux for all of the models, as shown in Figure 1 1 
for Region 2. Annual latent heat also decreased, but to a lesser extent. Snow pack evolution 
and runoff production showed little change. 

The SEM is a usefitl tool to interpret the behavior of the LSSs as a result of changes in the wind 
speed and shortwave radiation. To incorporate both the effect of shortwave radiation and wind 
speed, which is not used in the SEM, the SEM parameters were re-estimated using the 
experiment (2) LSS results. The total available energy forcing for each of the simplified models 
was adjusted to reflect the experiment (2) forcing conditions, and the SEM simulations were 
repeated. The change in energy partitioning resulting from the SEM calculation is shown in 
Figure 12 for Region 2. Despite the simplicity of the model and the lack of a representation of 
energy components, the model was able to capture the magnitude and direction of the changes 
in Bowen ratio between experiment (1) and experiment (2), on average. However, there is a 
tendency for the SEM to underestimate the magnitude of the decrease in latent heat flux 
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between experiment (1) and (2), with a subsequent over-estimation of the change in sensible 
heat. 

The coincidental change in both shortwave radiation and wind speed makes it difficult to 
interpret the results. However, changes in model performance during the winter months were 
due primarily to the change in the wind speed, because shortwave radiation in the Tome-Kalix 
basin is near-zero during the winter months. For Region 2, average wind speed increased 
between experiment (1) and (2), corresponding to a decrease of winter sensible heat for many 
models (Figure 13a). In Region 5, the average winter wind speed decreased by a similar 
magnitude, resulting in a consistent increase in winter sensible heat (decrease in magnitude) 

(Figure 1 3b). 

Figure 14 shows the change in the fitted parameters when estimated with the reduced 
shortwave radiation and changed wind speed model runs for Regions 2 (black) and 5 (gray). 
For both regions, the critical temperature gradient for switching between stable and unstable 
conditions became more negative for experiment (2). Consequently, the average effective 
resistance for Region 2 decreased, despite both positive and negative changes in the indivtdual 
resistances. The decreased resistance led to an enhancement of turbulent transfer of latent heat 
and a compensating increasingly negative sensible heat flux. For Region 5, both the stable and 
unstable resistances, and therefore, the effective resistance increased. This means that turbulent 
transfers were suppressed in experiment (2), leading to a decrease in the magnitude of the 
sensible heat flux, the opposite behavior to Region 2. Examination of the response of the other 
regions suggests that the Region 2 response is controlled by the change in wind speed and is 
little influenced by the presence or absence of overstory vegetation. 
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The average change in wind speed is near zero for Region 2 during the summer months, and 
modeled changes between experiment (1) and experiment (2) were primarily attributable to the 
decrease in shortwave radiation. Evaporation efficiency showed both positive and negative 
changes for Region 2 (Figure 14b). Because the evaporation efficiency is less than 0.5 for the 
majority of the models, most of the decrease in available energy is accounted for by a decrease 
in the sensible heat flux. In Region 5, the evaporation efficiency increased for nearly all models 
for experiment (2). The contrast with Region 2 suggests that the decrease in wind speed led to 
sensible heat consuming a smaller portion of available energy, and subsequently a smaller 

decrease in summer sensible heat. 

6. Discussion and Conclusions 

To explain the relationship among snow cover, turbulent exchanges and runoff production 
among the LSSs participating in the PILPS Phase 2(e) intercomparison, a simple equivalent 
model was devised to reproduce the seasonal and annual fluxes of 13 PILPS LSSs. This 
analysis shows that for the high-latitude Tome-Kalix River basin, the annual partitioning of 
precipitation and energy in the LSSs can be related to three primary parameters: snow albedo, 
effective aerodynamic resistance and evaporation efficiency. These parameters were estimated 
from LSS-simulated variables and used to parameterize the SEM for five regions within the 
Tome-Kalix River basin. On average, the SEM is better able to reproduce annual evaporation 
in areas without overstory vegetation (Regions 1 and 2). Because these areas are also the 
regions that receive the least annual radiation, this may be a reflection of a greater energy 
limitation on evaporation at these latitudes, in addition to vegetation controls. Over all models, 
the SEM is better able to represent the seasonal snow cycle in regions with overstory (Regions 
3, 4 and 5), as measured using the squared correlation coefficient. This may be due in part to a 
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greater sensitivity to the aerodynamic resistance in Regions 1 and 2. The estimated 
aerodynamic resistance parameters tended to be higher for Regions 1 and 2, and annual runoff 
production is most sensitive to the aerodynamic resistance in these regions. In Regions 3, 4 and 
5, annual runoff production is most influenced by both the effective aerodynamic resistance and 
the partitioning of energy into evaporation in the summer. Interannual variability in 
summertime evaporation is more influenced by interannual variation in model parameters in 
regions with overstory vegetation, then in regions without overstory. In all cases, however, 
interannual variation in available energy controls a larger proportion of the total variation m 

evaporation. 

The SEMs were somewhat successful in reproducing the observed LSS response to a decrease 
in shortwave radiation forcings and changing wind speed. Fitted parameters to the revised 
forcings indicate that a larger proportion of available energy, and thus a decrease in available 
energy, is taken from sensible heat. Changes in wind speed can counteract this effect, leading 

to subsequent decreases in latent heat. 
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Appendix A: Parameter estimation for the SEMs 

Using the method of Koster and Milly (1997), fta was found as the highest soil wetness for 
which the LSS produced monthly runoffless than 1 mm/day. The two parameters m R , and mar 
were then found as the average slope of total runoff (surface and subsurface runoff returned 
from the LSSs) versus LSS soil wetness for the data series above and below this threshold, 
respectively. A sample of the fitted regression lines are shown in Figure Al (a). 

The parameters for snow-free evaporation were found by partitioning the monthly time series 
of available energy and evaporation into snow covered and snow-free periods. When the snow- 
covered fraction was less than 10 percent the entire grid cell was considered snow-free. The 
parameters /)„ and E, were found by least squares regression of the LSS evaporation variables 
and available energy computed from the monthly net radiation, ground heat flux and 
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interception evaporation simulated by the LSS. A sample of the fitted regression lines are 
shown in Figure Al (b). 

The parameters of the snow model were derived using daily time series of the air temperature 
and surface pressure forcings and the surface temperature and sensible heat flux variables 
returned by each model for days on which the snow cover fraction as determined by each 
model exceeded 90%. Examination of plots of LSS simulated Qh versus p a ■ Cp air {T s -T air ), 
indicated that most models demonstrate two distinct modes of behavior, due to the utilization of 
corrections for the stability of the atmosphere, as shown in Figure A.l (c). In order to segregate 
these two modes of behavior, the unstable aerodynamic resistance, r a _ u „, was found as the slope 
of the least squares regression of Qh versus p m -CpjT t -Tj, for all days on which the bulk 
Richardson number was less the 0.0045. The critical temperature gradient, AT crit „ for 


switching between stable and unstable conditions was then found as. 


AT,, = 


_&( 98%) 



(Al) 


where Q h (98%) is the observed value of sensible heat flux that is exceeded on 98 percent of the 
days. The unstable aerodynamic resistance was found as the slope of the fitted line for all 
points for which T s - T air is greater than the critical threshold. 

Snow albedo, a, was calculated as the mean monthly value of snow albedo returned by each 
LSS, for each region. If a LSS did not return snow albedo, it was calculated as the mean land 
surface albedo returned by each model for months in which the snow cover fraction exceeded 
75 percent. The value was held constant at the June value for months in which the snow cover 
fraction did not exceed 75 percent. Similarly, a time series of snow covered area was 
determined as the monthly average value returned by each LSS. If snow cover fraction was not 
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returned by a LSS, the fraction was specified as 0 or 1 , depending on the presence or absence of 
snow. 


20/38 


Bowling et al. 2002. PILPS Phase 2(e) Part 3. Global and Planetary Change, to be submitted. 


Captions 


Figure 1 : Land cover map of the Tome-Kalix basin, showing the location of the five sub- 
regions used in the SEM analysis. 

Figure 2: The dependence of runoff and summer evaporative fraction on scaled soil moisture 
for all models for each of the five sub-areas: a) monthly evaporative fraction for summer 
months and b) monthly runoff. To clarify the relationship, the dots are connected in order o 
decreasing evaporative fraction (increasing runoff), not in temporal order. 

Figure 2: The dependence of runoff and summer evaporative fraction on scaled soil moisture 
for all models for each of the five sub-areas: a) monthly evaporative fraction for summer 
months and b) monthly runoff. 


Figure 3i Simple Equivalent Model Schematic 


Figure 4' Goodness of fit for the simple equivalent models versus the LSS simulations 
and b) annual evaporation bias, c) and d) annual soil wetness bias and e) and f) squared 
correlation coefficient of monthly SWE. The left hand column shows Regions 1 and 2, 
right hand column shows Regions 3, 4 and 5, coded by number. 


for a) 
and the 


Figure 5. Mean monthly SWE LSS simulated (solid line) versus simple equivalent model 
estimate (dashed line), for Region 2 (black) and Region 5 (gray). Letters correspond to the LSS 

listed in Table 2. 


Figure 6: Aerodynamic resistance parameters controlling turbulent fluxes from snow, for each 
region: a) stable case and b) unstable case. 


Figure 7: Average annual melt predicted by each bucket model as a function of fitted 
parameters: a) average aerodynamic resistance and b) average snow albedo. 


Figure 8: Average annual evaporation efficiency for no-snow condition, versus average annual 
available energy from the LSS: Region 2 (black) and Region 5 (gray). Arrows show magnitude 
of interannual variability (as represented by the standard deviation) for each quantity. 

Figure 9: Sensitivity of annual runoff production to bucket model parameters for Region 2 a) 
Amiual runoff with all parameters fixed, b) difference between original and fixed runoff, c) 
percent difference due to fixed R a , d) percent difference due to fixed evaporation efficiency, e) 
percent difference due to fixed albedo and f) percent difference due to fixed snow cover 

fraction. 


Figure 10: Sensitivity of annual runoff production to bucket model parameters for Region 5. a) 
Amiual runoff with all parameters fixed, b) difference between original and fixed runoff, c) 
percent difference due to fixed R a , d) percent difference due to fixed evaporation efficiency, e) 
percent difference due to fixed albedo and f) percent difference due to fixed snow cover 

fraction. 
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Figure 1 1 : Magnitude and direction of change in average annual sensible and latent heat flux 
from the original LSSs for Region 2 going between experiment (1) and experiment (2). 

Figure 12: Magnitude and direction of change in average annual sensible and latent heat flux 
from the SEM for Region 2 going between experiment (1) and experiment (2) forcings an 

parameters. 


Figure 13: Mean monthly change (experiment 1 - experiment 2) in incoming shortwave (gray 
dashed), sensible heat (solid black) and wind speed (black dashed) for the LSSs: a) for Region 
2 (no overstory) and b) for Region 5 ( overstory). 


Figure 14: Change in fitted parameters between experiment (1) and (2) for Region 2 (black) 
and Region 5 (gray): a) AT cr ,t, b) B ns , c) r a _ stable and d) r a unstable 

Figure Al : Sample fitting procedure for one LSS: a) Runoff, b) Evaporation and c) 
aerodynamic resistance. 
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Table 1 : Description of modeling sub-regions 


a auiv a . 

Group 



Annual 

Precipitation 

(mm) 

Percent 

Snow 

(%) 

Annual Air 
Temperature 

(°C) 

Mean 

Elevation 

(m) 

Dominant 

Vegetation 

1 

703 

48.4 

-2.1 

697 

Open/closed 
shrub land 

2 

1273 

55.1 

-2.2 

974 

Open shrubland/ 
bare ground 

3 

612 

43.1 

-1.4 

404 

Wooded 

grassland 

4 

746 

42.0 

0.0 

299 

Woodland 

5 

710 

38.5 

0.9 

126 

Wooded 

grassland 
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Table 2‘ Aver age simple equivalent model parameters for each LSS 

— Snow model Evap Runoff response 

Model A T cri , r a sl r a un a fL e q EE EE— 


A 

CLASS 

-2.5 

0.0013 

0.0134 

B 

ECMWF 

-2.2 

0.0034 

0.0516 

C 

HY-SSiB 

-4.1 

0.0020 

0.0060 

D 

ISBA 

-3.1 

0.0004 

0.0172 

E 

MECMWF 

-2.8 

0.0018 

0.0211 

F 

MOSES 

-2.5 

0.0057 

0.0436 

G 

MOSES-CEH 

-2.6 

0.0054 

0.0437 

H 

NOAH 

-2.9 

0.0106 

0.0251 

1 

SAST 

-2.1 

0.0056 

0.0367 

J 

SEWAB 

-9.2 

0.0001 

0.0079 

K 

SWAP 

-3.6 

0.0029 

0.0105 

L 

VIC 

-2.8 

0.0019 

0.0187 

M 

VISA 

-0.8 

0.0091 

0.0837 


0.72 

0.59 

0.67 

7.0 

0.4 

0.1 

0.75 

0.54 

0.79 

9.8 

0.5 

0.4 

0.62 

0.64 

0.26 

8.3 

0.8 

0.6 

0.80 

0.58 

0.64 

5.1 

0.4 

0.2 

0.83 

0.57 

0.78 

10.5 

0.8 

0.6 

0.81 

0.48 

0.94 

10.8 

0.9 

0.7 

0.81 

0.42 

0.91 

15.6 

0.8 

0.7 

0.44 

0.44 

0.48 

13.9 

0.9 

0.7 

0.68 

-0.06 

1.47 

18.8 

0.6 

0.5 

0.82 

0.48 

0.48 

9.6 

0.8 

0.7 

0.85 

0.64 

0.69 

8.4 

1.0 

0.5 

0.76 

0.53 

0.56 

8.5 

0.7 

0.5 

0.79 

0.48 

0.77 

13.3 

0.6 

0.4 


Mr2 

29.0 

40.8 

170.3 

102.3 

66.0 

10.8 

25.9 

13.9 
94.5 

159.3 
8.4 

152.8 

138.0 


Table 3: Parameter values 

used for sensitivity runs 




A T cr „ r a _sl r a _u n a 

Pns 

E 0 

Constant in time: 

-3.16° C 259 sm' 1 34.3 sm' 1 0.75 

0.49 

0.73 mm 

Monthly varying snow- 
covered fraction: 

0.99, 0.99, 0.99, 0.94, 0.60, 0.06, 0.00, 0.00, 

0.07, 0.56, 0.91, 0.96 
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Figure 1 : Land cover map of the Tome-Kalix basin, showing the location of the five sub-regions 
used in the simple model analysis. 
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Figure 2- The dependence of runoff and summer evaporative fraction on scaled soil moisture for 
alfmodels for each of the five sub-areas: a) monthly evaporative fraction for summer months and 
b) monthly runoff. To clarify the relationship, the dots are connected in order of decreasing 
evaporative fraction (increasing runoff), not in temporal order. 


26/38 


Runoff (mm/day) Runoff (mm/day) Runoff (mm/day) Runoff (mm/day) Runoff (mm/day) 



Bowling et al. 2002, PILPS Phase 2(e) Part 3, figures 



Figure 3: Bucket Model Schematic 
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Estimated Evaporation 




Estimated Evaporation 



Estimated Soil Wetness 



Cumulative probability 



Cumulative probability 


Figure 4- Goodness of fit for the simple bucket models versus the LSS simulations for a) and b) 
annual evaporation bias, c) and d) annual soil wetness bias and e) and f) squared correlation^ _ 
coefficient of monthly SWE. The left hand column shows Regions 1 and 2, and the right han 
column shows Regions 3, 4 and 5, coded by number. 
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Month 


Figure 5. Mean monthly SWE LSS simulated (solid line) versus bucket model estimate (dashed 
line), for Region 2 (black) and Region 5 (gray). Letters correspond to the LSS listed in Table 2. 
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Figure 6: Aerodynamic resistance parameters controlling turbulent fluxes from snow, for each 
region! a) stable case and b) unstable case. 


5 - 

a ) H iei LI E JD 

4 - 

WH 1 B LAX DEB J 

c 


o 

’5b 3 ~ 

M6BELA D J 

QJ 


CC 


2 - 

IBM BE LA J D 

1 - 

aar ec la d J 


b) 

OCI EH AD KC J 

M 

f-B D Hi C K J 


M Be in fO KC J 


B H« B<A DC 


B6 NIIEU AD C 

.Ml -l— 1 III TTTT -1—r* 


1 ■ ■ 1 1 m , h * 

10° 10 1 10 2 


R (unstable) 




Figure 7: Average annual melt predicted by each bucket model as a function of fitted parameters, 
a) average aerodynamic resistance and b) average snow albedo. 
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Figure 8- Average annual evaporation efficiency for no-snow condition, versus average annual 
available energy from the LSS: Region 2 (black) and Region 5 (gray). Arrows show magnitude 
of interannual variability (as represented by the standard deviation) for each quantity. 
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Figure 9: Sensitivity of annual runoff production to bucket model parameters for Region 1 a) 
Annual runoff with all parameters fixed, b) difference between original and fixed runoff, c) 
percent difference due to fixed R a , d) percent difference due to fixed evaporation efficiency, e) 
percent difference due to fixed albedo and f) percent difference due to fixed snow cover fraction. 
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Figure 10: Sensitivity of annual runoff production to bucket model parameters for Region 3. a) 
Annual runoff with all parameters fixed, b) difference between original and fixed runoff, c) 
percent difference due to fixed R a , d) percent difference due to fixed evaporation efficiency e) 
percent difference due to fixed albedo and f) percent difference due to fixed snow cover fraction. 
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Figure 1 1 : Magnitude and direction of change in average annual sensible and latent heat flux 
from the original LSSs for Region 2 going between experiment (1) and experiment (2). 
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Figure 12: Magnitude and direction of change in average annual sensible and latent heat flux front the SEM 
for Region 2 going between experiment ( 1 ) and experiment (2) forcings and parameters. 
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Figure 14: Change in fitted parameters between experiment (1) and (2) for Region 2 (black) and 
Region 5 (gray): a) AT cril , b) B ns? C) T a stable and d) I"a unstable 
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Figure Al : Sample fitting procedure for one LSS: a) Runoff, b) Evaporation and c) aerodynamic 
resistance. 






